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Abstract 

We consider the inverse problem of determining spatially heterogeneous 
absorption and diffusion coefficients D{x), from a single measurement 
of the absorbed energy S{x) = /a{x)u{x), where u satisfies the elliptic partial 
differential equation 

—V • {D{x)\/u{x)) + ia{x)u{x) =0 in O C 

This problem, which is central in quantitative photoacoustic tomography^ 
is in general ill-posed since it admits an infinite number of solution pairs. 
Using similar ideas as in [31], we show that when the coefficients ja,D 
are known to be piecewise constant functions, a unique solution can be 
obtained. For the numerical determination of /i, D, we suggest a variational 
method based on an Ambrosio-Tortorelli approximation of a Mumford- 
Shah-like functional, which we implement numerically and test on simulated 
two-dimensional data. 
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1 Quantitative photoacoustic tomography 

1.1 Introduction 

Photoacoustic tomography (PAT) is a medical imaging technique that combines 
electromagnetic excitation (in the visible spectrum) with ultrasound measure¬ 
ments. In a photoacoustic experiment, a translucent sample is illuminated by a 
laser pulse with wavelength A (which we assume to be fixed in this article). The 
absorbed optical energy leads to thermal expansion, generating an ultrasound 
pressure wave p(x, t) that can be measured outside the sample. 

Since the pulses used in PAT are very short, the complete energy is deposited 
almost instantaneously compared to travel times of acoustic waves, and the 
pressure wave p can be assumed to have been generated by an initial pressure 
1-L{x)^ that is, it satisfies the wave equation 

dttp{x,t) - c^{x)Ap{x,t) = 0 
dtp{x, 0) = 0 
p{x,0) = V-ix) 

and p\m^ where A4 denotes a measurement surfaee^ can be obtained from 
ultrasound measurements [13, 14, 28, 45]. 

By solving an inverse problem for the wave equation (see, e.g, [28] for a review 
of inversion techniques), these ultrasound measurements can be used to estimate 
the initial pressure 1-L{x). Since 

1-L{x) = V{x)£{x) = V{x)ii{x)u{x), (1.1) 

that is, 1-L{x) is proportional to the absorbed energy f(x), which is in turn 
proportional to the optieal absorption eoeffieient ja{x) at the applied wavelength 
A and the local fluence u{x) (the time-integrated laser power received at x), the 
initial pressure visualizes contrast in fi. The constant of proportionality r(x) is 
called Griineisen eoeffieient (or PA efficieney since it describes the efficiency of 
conversion from absorbed energy to acoustic signal) [13, 14]. 

In quantitative photoaeoustie tomography (qPAT), the goal is to apply PAT to 
determine (inhomogeneous) optical material properties of the sample (which are 
of diagnostic interest). To do so, an additional non-linear inverse problem for 
light transport has to be solved, since the fluence u{x) is inhomogeneous and 
itself dependent on the optical properties of the sample [13, 14]. 

While it is also possible to attack the acoustic and optical inverse problems 
simultaneously (see, e.g., [9, 24]), here, we assume that the acoustic part of the 
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problem has been solved successfully, i.e., that (possibly noisy) data ~ 1-L 
(with 6 signifying the noise level) are available. 

In this article, we utilize the diffusion approximation (which is valid in highly 
scattering media [13, 45]) of the radiative transfer equation to model the fluence 
distribution. It is, however, also possible to use a radiative transfer model for 
qPAT (see [17, 33, 38, 44, 46]), albeit at the cost of increased computational and 
analytical complexity. 

The diffusion model (in a Lipschitz domain Vl C M^) is given by 

—V • {D{x)\/u{x)) + /a{x)u{x) = 0 in {12) 

u{x)\dn = g{x). 

The parameter D is called diffusion coefficient. The Dirichlet boundary data 
g (which we assume to be continuous and known in this article) describes the 
illumination pattern. Note that this is a time-independent model, again due to 
the fact that energy is deposited almost instantaneously compared to time scales 
of the acoustic system. 

The difficulty of the inverse problem varies depending on which parameters of 
the model are assumed to be known. We present three inverse problems often 
considered in the literature. The hardest one is: 

Determine (/i,D,r) from measurements of . (P3) 

Bal and Ren showed in [6] that for arbitrary coefficients /i, D, P G IT^’^(^4), this 
problem is unsolvable, even if multiple measurements of H (with different known 
boundary illumination patterns g) are available. 

In [31], the authors showed that with a restriction to piecewise constant parame¬ 
ters^ unique reconstruction of all three unknown parameters /i, D, P from multiple 
measurements is possible (under a condition on the directions of Vu/., where 
Uk is the fluence of the kth illumination pattern). Furthermore, an analytical 
reconstruction procedure was suggested and implemented numerically, which, 
unfortunately, is relatively sensitive to noise. 

Alberti and Ammari [1] also established a unique reconstruction result, based 
on morphological component analysis (a sparsity approach), in a slightly more 
general setting (which assumes different degrees of smoothness of the coefficients 
and the fluence). They also provide numerical reconstructions, which were, 
however, not tested for noise sensitivity in the case of (P3). 

To simplify the problem, it is often assumed that the Griineisen coefficient T is 
known or constant, which implies that the absorbed energy can be estimated 
with ^ ^ £ (with S again denoting the noise level). It remains to solve 

Determine (ja^D) from measurements of £^. (P2) 

If only a single measurement of £^ is given, this inverse problem is also ill-posed 
(since it has infinitely many solutions pairs, see [12, 31, 39]). 
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However, in [1], the authors were able to recover /i (independently of the light 
transfer model used) from a single measurement of £ (again using a sparsity 
method^ assuming different degrees of smoothness of the coefficients and the 
fluence). 

In [6, 7] it was shown that this problem is uniquely solvable if two measurements 
of £ (corresponding to well-chosen boundary illuminations ^ 1 ,^ 2 ) are available. 
Numerically, this multi-illumination case was treated in [6, 22, 36, 39, 44, 48], 
using a multitude of different techniques, see also the review paper by Cox et al. 

[13]. 

The simplest case works under the assumption that the diffusion coefficient D is 
also known: 

Determine fi from measurements of 5^. (PI) 

This inverse problem has a unique solution even for a single measurement, which 
can be seen by substituting fiu = £^ in (1.2), providing the possibility to solve for 
u [8]. For other (numerical) approaches, cf. [13]. To the authors knowledge, this 
simplified problem is the only case for which practical viability (with experimental 
data) has been established both for phantoms [47] and biological samples [40, 41]. 

1.2 Contributions of this article 

In this article, we consider the problem (P2) using a single measurement for our 
reconstructions, a problem that is in general ill-posed. Similar to [31] (which 
treats (P3) with multiple measurements), a restriction to piecewise constant 
/i, D also proves to be useful for this problem. In fact, as shown in Section 2, 
when the parameters /i, D are piecewise constant functions (and noise-free data 
are given), the inverse problem (P2) can be solved uniquely, without any further 
assumptions. 

In Section 3, we present a variational model for the reconstruction of piecewise 
constant /i, D from noisy data £^ based on the Ambrosio-Tortorelli approximation 
of a Mumford-Shah-like functional. Compared to the two-step reconstruction 
process presented in [31] (which was introduced for (P3), but is applicable for 
(P2)), which first detects the regions where the parameters are constant, then 
reconstructs the parameter values from jumps of the data and its derivatives, 
this variational approach is much more robust with respect to noise. This is 
mainly due to the fact that the numerical approach presented in [31] requires 
almost perfect jump detection in the second derivatives of £^ to get reasonable 
estimates of /i, D (since the jumps have to form a full partition of the domain 
O), which is highly challenging in the presence of significant amounts of noise. 
This is not the case for the variational method presented here. 

Finally, a description of our implementation and numerical results can be found 
in Section 4. 
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2 Recovery of piecewise constant coefficients 

In this Section we show that piecewise constant parameters /i, D can be recovered 
uniquely from a single measurement of the absorbed energy 

In the following, let be a partition of 0 C into open sets and ji^D 

piecewise constant on that is, for open and (im^Djn G 

MM M 

m=l m=l m=l 

Furthermore, for /c G N, denote by 

Jk{f) =£t\ C I 5 is open and / G C^{B)} 

the discontinuities of a function / G L^{Q) and its derivatives up to k-th order. 
First, we show that coefficient discontinuities can be recovered from the data £. 

Proposition 1. Let ji^D he of the form (2.1) and £ = £{ia,D) satisfy (1.1),(1.2) 
weakly. Then we have 

Jo(/i)UJo(I)) = J 2 (f). (2.2) 

Proof. First, take an arbitrary open ball B with B D (Jo(/^) U Jo{D)) = 0 and 
notice that, by interior elliptic regularity, u G C^{B) and hence also £ G C^{B). 
It follows that 

Mti)UJo{D)DJ2{£). 

By the De Giorgi-Nash-Moser theorem [23, Theorem 8.22], we have u G (7^(0), 
so 

Jo(/i) = Jo(f). (2.3) 

Finally, as G (7^(r^rn), m = 1,..., M, (1.2) holds strongly in all and we 
get 

Dr/ji^u -|- = 0 in Llfri ^ .x 

= 0 in 

since jirn^Dm are scalars. Hence, A£ = ^ £ = in Um=i which shows 
that A£ cannot be continuous where ^ jumps (since u G (7^(f2)), hence 

Jo C J2{8). 

The Proposition follows from Jo(/i) U Jq = Jo(/a) U Jo{D). □ 

Remark 1. The proof of Proposition 1 also shows that jumps in fi and D have 
different effects on the data £. By (2.3), the jumps of fi can be recovered from 
those in £. Jumps of D (not coinciding with jumps of /i) on the other hand. 
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are smoothed in the data and have to be obtained from the derivatives of f. In 
particular, the discontinuities of Af suffice to find the jumps in D. 

Furthermore, note that from discontinuities of Vf (or, equivalently, |Vf p), one 
can, in general, only recover the part of the jumps of D where \/u{x) • u{x) ^ 0 
(where v denotes a normal vector to the hyper-surface of discontinuity in D), 
as may be continuous across parts of the jump set of D where \/u{x) is 
tangential. 

To see this, consider for instance the case of a partition consisting of a homoge¬ 
neous region Qi with (non-touching) smooth inclusions • • •, In this case, 
we have u G C^{Qrn) for all m = 1,..., M (due a result of Li and Nirenberg, 
see [29]) and the interface conditions 

■ jy = DiVu\q^ ■ jy, Vu\q^ ■ t = ■ t (2.5) 

hold point-wise on dflrn for m = 2,..., M, where z/, r denote vectors normal and 
tangential to dQrn (see, e.g., [31] for a derivation). Now, let 

M 

A := {x G dflrn I Vu{x) • iy{x) ^ 0}. 

m=2 

We have Jq{D) fl A C Ji{u) since, by the interface conditions (2.5), Vu cannot 
be continuous across Jo{D) D A (it changes length). On the other hand, if we 
take an open ball 5 C fl/c U (for some 2 <k < M) with B D Jo{D) D A = 0, 
we have, again by (2.5), 


( 2 . 6 ) 

since either = Di or = VizjoinB = 0- As we have u G C^(12/c), 

u G C^(12i), (2.6) implies u G C^{B) and therefore Ji{u) C Jo{D) D A, hence 

Ji{u) = Jo{D) n A. 

This shows that, in general, discontinuities in the second derivatives of the data 
£ have to be identified in order to get the whole jump set of D. 


Once the partition is known, the coefficients can be recovered 

from the jumps in f^ and the boundary values of u. 

Proposition 2. Let fi^D be of the form (2.1) with {Ltm)m=i known. Further¬ 
more, let £ satisfy (1.1),(1.2) for a known boundary illumination g. Then, ji 
and D can be determined uniquely from £. 

Proof Let x G d£lm H d£ln for some m,n £ {1, ..., M}. Since u G 

lilRy^x ^{y)\Qrn _ _ Tm ^2 ij\ 

lim;2^a,f (2:)|o^ /anU{x) /in* 
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Furthermore, take x G dQ D dClk (for some k G M}). We have 


lim 

y^x 


9{x) 


lik lim 

y^x 


u{y)\nk 

9{x) 


— l^k' 


Starting in and using (2.7) on all interfaces, we recover all iimi m = 1,..., M. 
Finally, to obtain m = 1,..., M, we use (2.4), that is. 


l^m 

Dm 


Ag(2/) 

^{y) 


for all y € Q,m- 


□ 

Remark 2. If g is not known. Proposition 2 can be used to determine the 
parameters /i, I) up to a constant. 


3 A Mumford-Shah-like functional for qPAT 

In Section 2, we showed that piecewise constant absorption and diffusion coeffi¬ 
cients /i, F), can be recovered from noise-free data £ by an analytical procedure 
which first determines the coefficient jumps, i.e., the partition and 

then the numerical values (/im^ In [31], such a two-step approach was 

implemented numerically (for the problem (P3)). 

In the presence of significant amounts of noise in the data, however, such an 
approach is infeasible, since it requires the detection of jumps of derivatives 
up to second order of the data . In particular, jumps may remain partially 
undetected (e.g., if the edge detection has to be restricted to first derivatives 
due to noise, see Remark 1 and the numerical examples in Section 4), leading 
to an incomplete estimated partition and therefore highly erroneous parameter 
estimates. 

To overcome this problem, we propose a variational approach favouring piecewise 
constant solutions that estimates the numerical values of piecewise constant /i, D 
and their jumps at the same time. 

There are multiple different methods for piecewise constant regularization of 
inverse problems. For instance, a popular class of methods is based on the level 
set method [32] or variations of it, e.g., [10, 11, 16, 42, 43]). This approach has 
been suggested for qPAT (using the radiative transfer model) in [17], however, it 
was only tested using multiple measurements. In this article, we use an Ambrosio- 
Tortorelli approximation of a Mumford-Shah-like functional (which was first 
suggested for electrical impedance tomography in [37]). The main advantage of 
this approximation is that we can utilize incomplete jump information (obtained, 
for instance, from jumps of the data or other means) to initialize the minimization 
procedure (see Section 3.3). Numerically, leads to faster convergence (and in 
some cases improved minimizers). Additionally, the number of segments does 
not have to be known in advance (in contrast to multiple level set methods) and 
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the minimization with respect to the jump indicator functions is a simple elliptic 
problem. 

We want to minimize the Mumford-Shah-like functional 
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+ v/ \Vii\^dx+^ [ \VDfdx 

^ Jn\K^ ^ Jq\Kd 

+ + I3dH^-\Kd), 


(3.1) 


where is the {N — l)-dimensional Hausdorff-measure and S{fi,D) the 

operator that maps the (unknown) parameters fi^D to the measurements £ 
satisfying (1.1),(1.2). The minimum is taken over all ji^D in suitable (that is, 
point-wise bounded from below and above) subsets of lT^’^(fl\i^^), 
and all closed sets Kd G (the jump sets of the coefficients /i, D). 

Functionals of this type, first introduced by Mumford and Shah for image 
denoising and segmentation in [30] , have been applied for a wide range of inverse 
problems, see, e.g., [20, 26, 27, 34, 35, 37]). The basic idea behind the functional 
is as follows. The discrepancy term 


£{fi,D)-£^ 
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forces minimizing coefficients /i, D to be close (in L^-sense) to a solution of the 
inverse problem f(/i,T)) = £^ (of which there are infinitely many), while the 
regularization terms 



|VD|2 dx + + PdH^-\Kd) 


force both the variation of /i, Z) (outside their jump sets Kd) as well as the 
hyper-surface area of Kd to he small. The regularization parameters 
control the amount of continuous variation the parameters may have, whereas 
control the complexity of the coefficient’s jump sets K^^Kd- 

Note that while in general, minimizers of this functional can have some continuous 
variation, they have to be close to piecewise constant for large values oia^^ao^ In 
fact, the Ambrosio-Tortorelli approximation of the functional (3.1), which will be 
introduced in Section 3.2, can also be used for piecewise constant regularization 
(in the limit oc, see Remark 5). 


3.1 Existence of minimizers 

Existence of minimizers of functionals like (3.1) (which are, in general, not 
unique) was first established in [18] for image segmentation and in [37] for 
regularization of non-linear operator equations. 










To ensure that £{ii^D) is well-defined and that T has a minimizer, point-wise 
bounds have to be enforced (see [25]), so we have to restrict the coefficients to 
= {f e : a < f <b a.e. in for some 0 < a, 6 < oo a priori. We 

also require the following Lemma, the proof follows the ideas in [19]. 

Lemma 1. The non-linear measurement operator 

£■- xl{£lf c L^{£lf ^ £{ii,D) =IJU{IJ,D), 

where u{ii^D) solves ( 1 . 2 ), is eontinuous. 

Proof. Let {fin^Dn) in 1/^(0)^. Furthermore, let Un : = 

?x(/in:^n) and u := u{ii^D) be solutions of (1.2) corresponding to the given 
coefficients. 

By the usual energy estimates for (1.2) (cf. [21, Chapter 6, Theorem 2]) we have 
||'^n|lwu2(o) ^ C(a, 6, fl) • 

Hence, there exists a (re-labeled) subsequence {uk)ken with Uk weakly in 
IF^’^(O) for some u G IF^’^(O). From the weak form of (1.2) we get for all 
V e VFi’^(fl) and A: G N 

{DV{uk-u),Vv)^2^^y + (/^K -^),^)l2(o) 

= {{D - Dk)VUk,Vv)^2^^y + ((/i-/i/c)^/c,u)^2(f^) 

— C'diD/c — D\\l2(^q;^ + ll/i/c — /^IIl2(0)) ll^llL2(aO) 11'^II W^oo (q) • 

Taking k ^ oo and using the fact that the left side acts as a bounded linear 
functional on Uk G IF^’^(O), we obtain for all v G 1F^’'^(0) 

{DX{u-u),Vv)^2t^Q)2 + (M(M-«),v)i2(n) = 0. 

By density of C W^’‘^{£l), this implies u = u. Since the same argument 

also holds for every subsequence of the original sequence {ujfjnen^ we get ^ u 
weakly in IF^’^(fl) and thus Un ^ u strongly in . 

Finally, since u G L^{£l) by the maximum principle, continuity of £ follows from 
ll/^n'^n “ /^'^IIl2(0) — ~ /^IIl2(Q) II'^IIl°°(0) II^’T' ~ ^IIl2(0) ll/^n || 1,00 ('q) . 


□ 

Remark 3. Since ||/||^2 < ||/|lLi(n) ll/llL°=(n) for all / e X^{Q,), the operator 
£: X^{Q)^ C L^{£l)^ L£‘{Q) is also continuous. 

Following [3, 37], we show the existence of minimizers of a weak form of X 
defined for coefficients /i, D in the space 5SV(fl), the special functions of hounded 
variation^ which may have jump discontinuities. For a quick introduction, see 
Appendix A. 
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We obtain the functional T : D SBV{^Y ^5 

[ \v^\^dx 

' JQ 


T{p,D) = ^\\£{p,D)-£^ 


L2(0) 


aD 

~Y 


[ |VD|" dx + H^-^[S{^i)] + Pd H^-^[S{D)], 

Jn 


(3.2) 


where S{f) denotes the approximate discontinuity set of a Lebesgue-measurable 
function / and Vji^VD the density of the absolutely continuous part (with 
respect to the Lebesgue measure) of the respective distributional gradients (see 
Appendix A). 

In this setting, the existence of minimizers can be established using the direct 
method and the SBV compactness theorem. 

Proposition 3. The weak Mumford-Shah-like functional T has at least one 
minimizer {ft, D) in fl SBV{fl)‘^. 


Proof Let (/i^,I)^) G n SBV{^})‘^ be a minimizing sequence of the 

functional (3.2). Clearly, we have for all n G N 

llMn|lz,oo(n) + _£ IV/XnP dx + <C<^ 

l|Dn|lL~(n) + iV/XnP dx + H^-^[S{D)] < C < (x^. 

By the tSSV-compactness theorem (see Section A), and using Lemma 1, first 
applied to then to the obtained subsequence of there exist {fi^D) G 
SBV(TlY ^ subsequence {lJik^Dk)ken (re-labeled) with 

{^lk, Dk) (A, D) strongly in 
(V/ife, V-Dfe) ^ (VA, V-D) weakly in 
H^-\S{p)] <liniinfif^-i[5(/Ufc)] 

k^oo 

H^-'^[SiD)] <liniinfiJ^-i[5(Z)fe)]. 

k^oo 

We also have (ft^D) G X^{Q)‘^ since is closed in Furthermore, 

it is a minimizer of P since the I/^(f2)-norm is weakly lower semi-continuous and 
E: X^(0)^ C L^(0)^ ^ LS{Tt) is continuous. □ 

Remark 4. If the minimizer {(i^D) has essentially closed jump sets S{fi)^S{D) 
(that is, taking the closure of the jump sets does not increase their -measure) 

the quadruple (^fi^ S{fi)^ S{D)^ also minimizes the strong Mumford-Shah 

functional (3.1) among all /r G X^(f]) H D G X^(f]) 

and K^^Kd closed. While essential closedness of jump sets of Mumford-Shah- 
minimizers can be established under certain conditions on the discrepancy term 
[18, 25, 37], this is out of the scope of this article. 
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3.2 Approximation 

Using the approach of Ambrosio and Tortorelli (cf. [4]), one can obtain functionals 

: xl{nf n X n ^ m 


that approximate the functional X from (3.2) and are easier to minimize: 


Xe{lJ.,D,Vf„VD) = 1: £{ld,D)-S 


L2(n) 


+ ^ / ivl + C^)\Xii\^dx+^ f {v% + CD)\XDfdx 

t/ O J ^ 

In ^ dx 

I^D J ^ dx. 


(3.3) 


It is well-known that, for all a, /3 > 0 and ( = o(e), 

aJ^{v^+0\Xf\^+l3 (^e|Vu|2 + dx dx+l3H^-^[S{f)] 


in the sense of F-convergence in (formally, the latter has to be extended 

to an additional variable, see [4, 37]). Since the data term is continuous in 
L^{Q) and F-convergence is stable under continuous perturbations, we thus 
get F-lime^ 0 ‘^e = and therefore l/^(0)-convergence of (a subsequence of) 
J^£-minimizers (which can be shown to lie in a compact set [4]) to J^-minimizers. 

The existence of a minimizer of (3.3) can be established using the direct method. 
Proposition 4. The functional has at least one minimizer {ft^ D^v^^vd) in 


Proof. Fg is coercive with respect to the semi-norm 

|F|vvu2(o) + hnWw^’^iQ) + W'^DWw^’^in) ^ 

which, combined with the restriction (/i, F) G X^(Q)^, shows that a minimizing 
sequence is bounded in IT^’^(O)^, so it contains a subsequence that converges 
weakly to {ft, D^v^^vd) in kF^’^(O)^ (and thus strongly in 1/^(0)^). Moreover, 
note that the spaces X^{ft) and Xo(^4) are closed in L^(0). The Proposition 
follows since the discrepancy term is continuous in and the regularization 

terms are weakly lower semi-continuous in kF^’^(Q)^ (cf. [15, Theorem 1.13] for 
the Qf-terms). □ 

Remark 5. If we vary a^, an with e and take 




e^O 


oc 
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the minimizers of converge to piecewise constant (in 5SV-sense, that is, with 
gradient vanishing almost everywhere) functions (/i, D) that minimize 










among all piecewise constant A proof (for the single variable case) can be 

found in [4]. This explains, in light of Section 2, why this type of regularization 
is useful for qPAT. 


3.3 Minimization 

In this subsection, we proceed formally, i.e., without proving convergence, ex¬ 
istence of minimizers of sub-problems and that the necessary derivatives and 
adjoints exist. For the minimization of (3.3), we suggest the following alternating 
directions approach: 

Algorithm 1. 

1. Choose parameters Pd >0. 

2. Find an initial estimate K of the edge set Jo(/i) U Jo{D) using edge 
detection. 

3. Initialize with fjP = = 1 and = 1 — . 

4. Iterate until convergence: 

(i) = argmin^ D,v^,vD 

(ii) = argmin„^ 

We now explain the steps in more detail. 

Initialization using edge detection 

In our numerical simulations, it became clear that using Proposition 1 to obtain 
a reasonable estimate K of Jo(/i) U Jo{D) leads to faster convergence and better 
minimizers than simply taking iC = 0. 

Given the noisy measurements , we have to estimate the discontinuities (or 
edges) of \\/ £^ p and \A£^ \. Since the jumps in these functions are of 
multiplicative nature (proportional to the local value of rt or • z^), it is 
advantageous to apply a logarithmic transformation prior to edge detection (to 
obtain constant contrast), see [31]. 

Due to noise, the data has to be smoothed prior to taking derivatives (turning 
discontinuities into areas with large gradients). We smooth by convolution 
with a Gaussian, since this approach has the advantage that differentiation and 
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smoothing can be done in one step (by differentiating the low-pass filter instead 
of the function). We have 


Gcr(x) = ( 27 r) 2 cr ^ exp 


—X 


(VG<,)(x) = (27r)-Tcr-(JV+2)(_2.)gxp 

|2 


(z^ 

\2a^ 


(AG<,)(rc) = (27r)-^(7-^+2 exp 

Similar to [31], we proceed as follows: 

1. Detect edges Kq of /o := log * 1b Go-q) \neB 

2 . Detect edges of /i := log 




{£^*1b VG<,,)pV7))l 


(n\Ko)eB 


3. Detect edges K 2 of /g := log ( (s^ * 1 b AG^^) ) l(n\(KouKi))eB 

4. Take K = K 0 UK 1 UK 2 


Here, 1b acts as a cut-off function for the filters (e.g., using B = the p-norm 
ball of radius p). The operation AO B = {z ^ Vi \ {B z) <Z A} denotes set 
erosion, which is performed to avoid multiple detection of edges (by ensuring 
that the cut-off filter does not intersect with already detected edges or the 
outside of the domain). Note, however, that for large B this may lead to parts 
of edges (close to the domain boundary or already detected edges) not being 
detected. The scalar 7 is a minimal value enforced for | V f ^ p (to avoid creating 
singularities at zeros). 

The edge detection itself is performed by applying thresholds ^ 05 ^ 15^2 to the 
functions |V/o|, |V/i|, IV/ 2 I (taking the super-level sets as edge sets). Note 
that in contrast to [31], no complete segmentation is necessary and cruder edge 
estimates suffice, that is, the detected edges don’t have to be reduced to thin 
curves. 


Iteration 

In step (i), to find, for fixed v^^vd G ^ IT^’^(ll), 

argmin Be{p,D) 

we use Gauss-Newton-minimization. That is, in every iteration of an inner loop, 
we linearly approximate S{pk ^ Dk ^ sd) ^ £{pk,Dk) ^ {pk, Dk){s^, sd) 
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and take update steps 


s = 


arg mm 
s^,,SDew^^^in) 


^k) “ 1 “ ^ if^k’) Dk){^fi') ^d) ^ 


L^(n) 


+ ^ / ('^^ + C/x)|V(/i/c + dx + ^ f {v% ^ (D)\^{Dk sd)\‘^ dx 
j vt J vt 


/n ^ Jn 

which we then project into so we update with 

(/^/c+i? -^/c+i) — {f^ki -^/c) “h (5 A n) V 6 . 


(3.4) 


A straightforward calculation shows that a minimizer s = (s^^sd) of (3.4) 
satisfies the weak form of 

DkT Dk){s^, sd) + 5 ^, sd^ 

= —£ [f (/i/c, D/g) —£ ] — /i/c, 

5 d^Dk) 

with homogeneous Neumann boundary conditions for s^^sd and 

^ Vk“^’^(f^), X -V • (crVx). 

The linear operator £'{jii^D) and its (formal) adjoint £'{ii^DY are given by 


(3.5) 


£'{fi,D){s^,SD) = ^ • (sdXu)] 

S'ifi, D)* t={tu- uL;^^{fit),-Vu ■ VL-\ 

where u = i4(/i, D) and 


(3.6) 


Lf,,D ■■ Wo’^in) W-^’^{n), x^-V- {DVx) + iix. 

By introducing auxiliary variables ^ 1 ,^ 2 , ^3 (and taking fi = D = Dk): the 
equations (3.5),(3.6) can be re-written as the system of (weak) PDE 


{Sfj^u + iiyi + iiy2)u - uys + s^ = R^, 


ds^dn 

dv 


= 0 


(3.7) 


—Vu • V ^3 + Sd = Rd, —^ ^ 

L^i^dV?. = y{s^u^ yyi +/i^2), ^slao = 0 
Lfi,D ^2 = V • (sdXu), y2\dn = 0 
L^^dVi = -s^u, yi\sn = 0 

which is amenable to discretization as a sparse matrix using, e.g., the finite 
element method (FEM). 
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In step (ii) of the outer loop we have to find, for fixed jJ,, D £ -^a(^) H W^’‘^(yi), 


argmin F^(v^,vd). 

v^,VDexl(Q.)nw^’'^{Vt) 


The minimizer {v^,vu) satisfies the linear equations 

^-2/3eA +(a^lV/rp + ^)Id^ = 0 

(^-2/3eA + (a,,|VD|2 + ^)Idj = = 0 


which can be solved separately and implemented numerically using FEM. 


4 Implementation and numerical results 

In this section, the proposed algorithm is tested on two sets of simulated two- 
dimensional data (with different parameter range and level of detail). We also 
vary the noise on the data since the reconstruction quality strongly depends on 
the noise level. 

The data (see Figures 1 and 2) was generated in the diffusion model (1.2) using 
self-written (linear-basis) finite element code in MATLAB. For both examples, 
we took Vt = [0, 5]^ and used a uniform boundary condition g = 1. The simulated 
data were generated on a (400 x 400)-grid and then down-sampled (by averaging) 
to (200 X 200) to avoid inverse crime. After that, Gaussian noise with different 
intensities (standard deviations of 0.5% and 10% of the average signal value 
f (x) dx) was added to the data. 

The edge detector described in Section 3.3 was implemented by finite difference 
approximations of |V/o|, |V/i|, IV/2I using central differences inside and one¬ 
sided differences near the boundary of the domain. The functions /i, / 2 , fs were 
calculated by convolution filtering with the MATLAB function imfilter. For all 
examples, we used a square cutoff function, i.e., B = B^ with p = 2\(7k\ (where 
k = 0 , 1 , 2 ). 

The edge detector is used to detect jumps in the derivatives of the data 
up to second order (to obtain an initial estimate of the parameter jump set 
Jo(/i) U Jq{D)). Since this process is highly sensitive with respect to noise, we 
varied the edge detection procedure subject to the amount of noise in the data. In 
the noise-free examples, we estimated the jumps of all three functions /o, /i, /s, 
that is, jumps of derivatives of up to second order. We restricted the jump 
estimation to /o, /i for the low-noise examples (i.e., jumps of derivatives up to 
first order) and /o in the high-noise examples (only jumps in the data £^ itself). 

To obtain the parameters /i, D given an initial estimate of their combined edge 
set Jo{p) U Jo{D)^ we used the Ambrosio-Tortorelli approximation (3.3) of the 
Mumford-Shah functional introduced in Section 3. To minimize the functional. 
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(c) S{ii,D) (d) S{ii,dY (10% noise) 

Fig. 1 : Example A. Parameters /jl^D and FEM-generated data 


we used Algorithm 1, iterating, for all examples, the outer alternating-directions 
loop until the functional changes by less than 0.01% and the inner Gauss- 
Newton loop until the functional value changes less than 1%. We directly 
implemented the systems (3.7) and (3.8) using (self-written) linear-basis finite 
element code. The implementation is not fully conforming, that is, where 
necessary we used conversions between piecewise linear and piecewise constant 
functions for simplicity. The discrete systems were solved using the standard 
MATLAB sparse equation solver mldivide. 

For all examples, we chose a = 0.01, 6 = 3 and e = 0.01. The other reconstruction 
parameters used (which were selected by hand) are listed in Tables 1 and 2. 

Reconstruction results and error profiles at different noise levels can be seen 
in Figures 3 and 4. In both examples, the noise-free reconstructions are very 
accurate and contain mostly smoothing error. In the low-noise reconstructions, 
due to the fact that more regularization is necessary, some of the parameter 
variation is underestimated. In the high-noise examples, most detail in D is lost 
since a lot of regularization is required to get reasonable results. The fine detail 
in /i can, however, still be recovered very accurately in both examples. 
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(a) n (b) D 





0.5 

0.45 

0.4 

0.35 

0.3 

0.25 

0.2 

0.15 

0.1 

0.05 

0 




(c) S{/j,,D) (d) (10% noise) 

Fig. 2 : Example B. Parameters and FEM-generated data £{jii,D). 


Remark 6. In some of the examples, minimization of (3.1) using log-parameters, 
i.e., the mapping £: (log/i, log 1)) f(/i,E)) instead of f, gave slightly better 

results. Note that this approach leads to a different linearization and therefore 
also a different system (3.7) to be solved in every step. However, to keep the 
presentation simple, we decided to use the functional as presented in Section 3 
for our numerical experiments. 

We also tried to incorporate the Griineisen coefficient E as an additional unknown 
into the reconstruction process, that is, to solve the problem (P3). Unfortunately, 
this proved to be highly unstable, even if the initial edge set was detected perfectly. 
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^/i, 

OLD 

/ 3 m 


C/x 

Cn 

Ex. A (0% noise) 

10-2 

10-^ 

10-5 

10-® 

10-5 

10-4 

Ex. A (0.1% noise) 

10-2 

10-5 

10-5 

10-3 

10-5 

10-3 

Ex. A (10% noise) 

1 

5•10-3 

10-5 

5•10-5 

10-5 

10-3 

Ex. B (0% noise) 

10-2 

10-^ 

10-5 

10-3 

10-5 

10-5 

Ex. B (0.1% noise) 

10-2 

10-5 

10-5 

10-3 

10-5 

10-3 

Ex. B (10% noise) 

1 

10-3 

10-5 

10-^ 

10-5 

10-3 


Table 1: Parameters used for Figures 3 and 4 (for the functional (3.3)). 



0-0 

0-1 

0-2 



^2 

7 

Ex. A (0% noise) 

0.5 

0.5 

0.5 

0.1 

0.1 

0.1 

10-4 

Ex. A (0.1% noise) 

0.5 

2 


0.05 

0.06 


10-5 

Ex. A (10% noise) 

1.5 



0.05 




Ex. B (0% noise) 

0.5 

0.5 

0.5 

0.1 

0.1 

0.1 

10-5 

Ex. B (0.1% noise) 

0.5 

2.4 


0.05 

0.06 


10-^ 

Ex. B (10% noise) 

1.5 



0.03 





Table 2: Parameters used for Figures 3 and 4 (for edge detection). 


18 





















Fig. 3: Reconstructions A. Columns: 0%, 0.1% and 10% noise. Rows: Estimated 
edge set, reconstructed parameters /igst^ ^est s^nd error profiles for /iest^^est- 
The error profiles show values of the true and reconstructed parameters along the 
image diagonals. The color axis in the images of the reconstructed parameters 
was fixed to the same range as the true parameters shown in Figure 1. 
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Fig. 4: Reconstructions B. Columns: 0%, 0.1% and 10% noise. Rows: Estimated 
edge set, reconstructed parameters /iest, ^est and error profiles for /iest^^est- 
The error profiles show valued of the true and reconstructed parameters along the 
image diagonals. The color axis in the images of the reconstructed parameters 
was fixed to the same range as the true parameters shown as in Figure 2. 
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A Special functions of bounded variation and 
the SBV-compactness theorem 


This section briefly introduces the notion of 5SV-functions and their compactness 
theorem. For a more comprehensive presentation with proofs, see, e.g., [5]. 

For a function / G 1/^(0) with distributional gradient Df^ we deflne its total 
variation by 

|D/|n := sup{{Df,ifi) \ ip e ||¥’||ioo(n;RiV) < 1} 

The space SV(0), consisting of all l/^(0)-functions with flnite total variation 
(i.e., of bounded variation), is a Banach space with the norm 

ll/ll23V(Q) •“ \\f\\L^(n) + \Df\n- 


Note that by the Riesz-Markov representation theorem, functions / G L^{Q) are 
of bounded variation if and only if Df is a finite vector Radon measure. The 
measure Df can be decomposed into three parts 

Df = D-f^D^f^D^f. 


D^f is the part of Df that is absolutely continuous with respect to the Lebesgue 
measure i.e., 

D^f = VfC^ 

for some integrable density function V/. The jump part D^ f is concentrated on 
the jump (or approximate discontinuity) set S{f) deflned by 

S{f) := {a; e O I /“(a;) < /+(a;)}, 


where, denoting with Bp{x) the ball centered at x with radius p. 


f+{x) := inf 


/ (a:) := sup 


t G M I lim ■ 

p^O 

t G M I lim 


C^{{yeB,{x)\f{y)>t}) 

C^{{y e Bp(x) I fiy) < t}) 
CN{B,{x)) 



Furthermore, for ^-almost all x G £'(/), there exists a unit normal vector 
u{x) and we have 

The remaining Cantor part D^f is concentrated on a subset of \ S{f) with 
intermediate Hausdorff dimension between N — 1 and N. 

A function / G is a special function of bounded variation (or / G 5SV(f])) 

if D'^f = 0, that is. 


Df = \7fC^ + (/+ - 


Furthermore, we have the following compactness theorem due to Ambrosio [2]: 
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Theorem 1. Let fn he a sequenee in SBV{Ll) with 


ll/nllLoo(n) + |V/„r dx + H^-^[SiU)] <M<^ 

for all n gN and some p > 1, M > 1. Then there exist a suhsequenee {fnk)kef^ 
and f G SBV{ft) with 

fuk -> / strongly in 

V/ weakly in LP{Q.-R^) 

< liminf 

k^oo 
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